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Abstract 

We present evidence for a novel finite-temperature phase transition in the 
phason elasticity of quasicrystals. A tiling model for energetically stabi- 
lized decagonal quasicrystals has been studied in an extensive series of Monte 
Carlo simulation. Hamiltonian (energetics) of the model is given by nearest- 
neighbor Penrose-like matching rules between three dimensional unit cells. A 
new order parameter and diagnostics have been introduced. We show that 
a transition from locked phason to unlocked phason dynamics occurs at fi- 
nite temperature. In the unlocked phase, phasons can be thermodynamically 
excited even though the quasicrystal is energetically stabilized at low temper- 
atures. 
64.70.Kb 
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I. INTRODUCTION 



Quasicrystals exhibit two distinct types of low-energy elastic ( hydrodynamic ) exci- 
tations — phonons and phasons§. In this paper we consider the temperature-dependent 
behavior of the phason elasticity in three-dimensional quasicrystals with decagonal sym- 
metry: a solid that can be described as a stack of periodically spaced planes which each 
exhibit ten-fold symmetryil. We present evidence for a novel finite-temperature phase tran- 
sition in the phason elasticity, apparently analogous to the pinning transition found in the 
one-dimensional Frenkel-Kontorova model (FK model)i. 

Quasicrystals are new types of solids which have a discrete point group symmetry that 
is forbidden for crystals such as five-fold symmetry in two- dimensions and icosahedral sym- 
metry in three-dimensionsS. These quasicrystals possess a long range translational order 
known as quasiperiodicity. The recent experiments on AlCoCu showed the existence of 
thermodynamically stable decagonal quasicrystals!. What makes the quasicrystals stable? 
Two possibilities that have been debated are energetic stability and entropic stability!. In 
energetically stabilized quasicrystal model, microscopic interaction energy has its minimum 
when atoms are arranged in a quasiperiodic structure. Such interactions ensure that the 
low-temperature equilibrium state is quasicrystalline. In entropically stabilized quasicrystal 
model, the entropy that arises from thermodynamically excited atomic relocations (specifi- 
cally, phasons) makes the quasicrystal stable. In this model, quasicrystals are not stable at 
low temperature. 

In this paper, we consider the phason dynamics of energetically stabilized tiling model 
for the decagonal phase. The interactions are prescribed so that quasicrystals remain stable 
as the temperature T approaches zero. A tile is an idealization of presenting a cluster of 
atoms in a real material. The energetics are mimicked by nearest neighbor matching rules 
which are generalizations of the Penrose edge-matching rules for two-dimensional tilings. 
We assign a finite energy to each mismatch in a given configuration. This guarantees that 
the state of lowest energy is a perfect quasicrystal. This is to be contrasted with a random 
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tiling model, used to represent the limiting case of entropically stabilized quasicrystals, in 
which the same energy is assigned to every configuration and a quasicrystal symmetry is 
favored due to the high entropy. 

Phonons and phasons are low energy hydrodynamic modes associated with quasiperi- 
odic broken translational symmetry. At long wavelengths, phonons correspond to uniform 
translations, and phasons correspond to rearrangements of atoms from one perfect quasicrys- 
tal lattice to another. In a tiling picture of quasicrystals, the phason degrees of freedom 
correspond to the rearrangements of tiles. Finite wavelength phasons produce rearrange- 
ments which violate the matching rules and hence cost finite energy. If w(x) is the phason 
field, then fixed phason strains produce a number of mismatches proportional to |Aw|. 
Consequently, the elastic energy is F ~ |Aw|, a nonanalytic form. We shall call a finite 
temperature state in which the elastic energy has this form a "locked phase"!, since the 
phasons cannot be thermodynamically excited (no phason Debye- Waller contribution). An 
alternative type of quasicrystal is described by the continuum density wave picture in which 
the elastic free energy grows F ~ (Aw) 2 . This density wave picture corresponds to a distinct 
elastic phase which we shall refer to as the "unlocked phase". In this phase, phasons have 
thermal excitations analogous to phonons. Historically, the locked phase has been assumed 
to be characteristic of any energetically stabilized quasicrystal, whereas entropically stabi- 
lized quasicrystals are, by definition, in the unlocked phase with square gradient elasticity. 
Our point here is to show that this view is not correct in general. Rather it is possible in 
energetically stabilized quasicrystals to have a novel elastic phase transition from a locked 
phase at low temperature to unlocked phase at high temperature as had been speculated by 
Socolar et. a/.i. We note that this behavior is different from 2D where both energetic and 
entropic model have unlocked elasticity at any finite temperaturei'0. Hence, the observa- 
tion of unlocked elastic behavior (phason Debye Waller effect ) at finite temperature is not 
a proof of entropic stability. 

The organization of the paper is as follows. In sec. II, we define our model for studying 
decagonal phase quasicrystals. In sec III, we present the characteristics of the system we use 
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for Monte Carlo simulation. The algorithms of Monte Carlo runs are also presented. Sec. 
IV presents diagnostics for studying the phase transition in phason elasticity. We introduce 
the notion of "lane width" and "trail magnetization" to check if the system is in the locked 
phase. The results of the Monte Carlo simulation are presented here. Sec. V, the conclusion, 
summarizes our results and discusses the potential connection to the experiment. 

II. DEFINING MODELS 

Our model for 3D quasicrystals with decagonal symmetry has two types of unit cells: 
skinny and fat. Each cell is the shape of a prism with rhombic cross-section whose upper 
and lower faces are the shape of Penrose rhombuses. In the zero temperature ground state, 
each layer viewed along the ten-fold axis resembles a perfect Penrose tiling. The 2D Penrose 
edge rules are replaced by "side face" rules, referring to tile faces that join tiles in the same 
layer. Side-faces of the unit cells have single or double arrow marks as shown in Fig. 1 so 
that the tiling consistent with matching rules — joining side-face to side-face according to 
the same arrow marks, is the perfect, quasiperiodic Penrose tiling. 

We assign mismatch energy e\ for a violation of single arrow marks. Double arrow 
mismatches are not allowed (mismatch energy for double arrow mismatch is infinity) in 
order to study purely phason dynamics(see discussion below Eq. (3)). The 3D decagonal 
model is introduced by considering a stack of these layers. We fix the interlayer interaction 
energy to have its minimum when configurations of two layers are identical (fat tile directly 
on fat tile, skinny directly on top of skinny). We assign stacking direction face mismatch 
energy e 2 when the vertices of an upper [lower] face of a unit cell do not coincide with those 
of a lowerfupper] face of an upper[lower] unit cell. The energy of a tiling is defined as a sum 
of intralayer interaction energy (arrow mismatch energy) and interlayer interaction energy 
(stacking direction face mismatch energy): 

^ ~ -^intralayer ^ ^interlayer 
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where 



£ intra 



a2 £ intra + 51 ^inter (-0 

side faces upper, lower faces 



for an arrow matched side face 

£1 for a single arrow mismatched side face 

oo for a double arrow mismatched side face, 



^inter 



and 

for an upper [lower] face which coincides with 

lower [upper] face of an upper [lower] unit 
€2 for an upper [lower] face which does not coincide 
with lower [upper] face of an upper [lower] unit. 

Henceforth we will refer to intralayer mismatches as "arrow mismatches" and interlayer 
mismatches as "face mismatches" . 

In the finite temperature Monte Carlo, we introduce " flips " of single hexagonal prisms 
within a layer. There are two kinds of hexagonal prisms which we shall call D-type hexagonal 
prisms and Q-type hexagonal prisms, respectively. A D[Q]-type hexagonal prism consists 
of two fat cells and one skinny cell [one fat cell and two skinny cells] and its top view is a 
D[Q]-type hexag on0. Each hexagonal prism can be uniquely identified by the vertex at the 
center of the upper hexagon. Fig. 2 shows hexagonal prism before and after "flips". Solid 
circles indicate the center vertices of the upper hexagons. 

A hexagonal prism flip from a ground state makes two arrow mismatches and 6 face 
mismatches. Hence, the net energy cost is 2e± + 6e2- 

For most of our MC runs, we impose a "stacking constraint": For a given hexagonal 
prism, flipping is allowed only when there is a hexagonal prism in the adjacent upper [or 
lower] layer whose lower [upper] hexagon has a boundary that coincides (ignoring center 
vertex) with the boundary of upper [lower] hexagon of the given hexagonal prism. This 
constraint is necessary for "interlayer flips" 
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Let us define phason variables in this model. The position of any vertex of the tiling 
with fixed origin (by choosing a vertex as the origin) is expressed as 

5 

x = n a el, (2) 

where e|, = (cos ^p, sin ^p, 0) for a = 0, 1, • • • , 4 and (0, 0, 1) for a = 5, n a is the number 
of steps in direction on a continuous path to the vertex at x from the origin along the 
edges, counting negatively when going in the — e|. We define a complementary basis e^, 
such that the vectors e a = © e][ are linearly independent in 6D hyperspace. Then, for 
each vertex x whose position is given by Eq. (2), we define a perp-space position vector x -1 
by 
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x± = E n « e a ( 3 ) 

where = (cos ^p, sin ^p, 1) for a = 0, 1, • ■ ■ , 4 and (0, 0, 0) for a = 5. These vectors span 
a 3D space called 'perp-space' which can be further decomposed as the product of a 2D- 
'phason space' (which is incommensurately oriented with respect to the hyperspace lattice) 
and an lD-'discrete space' (which is commensurately oriented). In our model in which double 
arrow mismatches are disallowed, the discrete space component of x -1 is restricted to four 
consecutive integers. Phason variables w are to be defined as a smoothed function w(x), 
constructed as an average of the phason space projection xi of the perp-space position 
vector X- 1 , in some neighborhood of radius a 3> 1. 



III. MONTE CARLO SIMULATION 

We use a thermal Monte Carlo methods based on Metropolis importance sampling 
scheme^. The basic Monte Carlo move is: 1) Randomly select a hexagonal prism. 2) if it 
satisfies the stacking direction constraint, flip it according to a probability; p = exp(—(3AE) 
if AE > 0, p = 1 otherwise, where AE is the energy cost for performing the flip. We 
used an approximant to a 2D perfect Penrose tiling for a layer which has the smallest phason 
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strain at a given system size with periodic boundary conditions. To get a periodic tiling, 
we use a periodic uniform dual grid which is the intersections of 4D lattice planes of the 5D 
hypercubic lattice with a 2D hypersurface spanned by two vectors, 

C3 l = (P, 0, -P, -Q, Q), u 2 = (0, P, Q, -Q, -P), (4) 

where P and Q are integers. A sequence of approximants are obtained by taking P = 
Pfc, Q = Fk-i, where is kth Fibonacci number (P = 0, Pi = 1). Then the basis vectors 
of the periodic tiling are given by 

L (1) =r k {*LZ1, -isLa^), L< a > =r* (0, 2sin^), (5) 
and the uniform phason strain E to get the kth periodic approximant is 

I 1 0^ 



E = -1 ft+i T 



k+l-2k 



V 2sinf tj 



(6) 



When we calculate perp-space position vectors, to get rid of the effects due to the uniform 
phason strain E, we replace in Eq. (3) by — Ee^ . The number of tiles in a layer is 
given by N k = F 2 k+i + 2P 2 fc . Since our main interest is the properties of 3D objects, the 
number of layers L z of the systems is proportional to the size of a layer (L z = P& for the 
jfeth approximant). Systems of 228(iV fe x L z = 76 x 3), 995(199 x 5), 4168(521 x 8) and 
17732(1364 x 13) tiles are used. We also study effect of varying L z for fixed L 1 and L 2 . For 
most of the models presented here, we have assigned interaction energy strength e\ = 1 and 
C2 = 1/3 in Eq. (1), corresponding to equal interlayer and intralayer energy cost for flipping 
an isolated hexagon beginning from the ground state. 

We start out with an ordered configuration, with the minimum matching rule violations 
necessary to construct the periodic approximants. Data are taken following a heating se- 
quence. We studied the time evolution of the data to check whether the system has been 
equilibrated. To ensure statistical independence, we measured the temporal correlations 
and most case measurements are separated more than correlation time0. Some cooling runs 
have been performed from above T c to below T c to check the presence of hysteresis effect. 



The data so obtained agree in high accuracy with those from heating studies. As a test 
of the algorithm, we assigned e 2 = without the stacking direction constraint and get the 
known 2D Penrose mo delS results. 



IV. DIAGNOSTICS FOR PHASON ELASTICITY TRANSITION AND RESULTS 

OF MONTE CARLO SIMULATION 

In the unlocked phase, where the continuum density wave picture is applied, the general 
form of the phason elastic free energy (up to quadratic in <9w) is 

F = 2 I ddxK ijkidiWjd k wi, (7) 

where Wj is the jth component of phason variable and is an elastic constant tensor. 
For our model, due to the decagonal symmetry, the free energy of Eq. (7) reduces to the 
form 

F = §/ d 3 x K^w,) 2 + (d yWl ) 2 + (d x w 2 ) 2 + (d y w 2 ) 2 } 

+K 2 [(d z w l ) 2 + (d z w 2 ) 2 ] + K 3 [d xWl d y w 2 - d yWl d x w 2 ) (8) 

where Ki,K 2 and K% are elastic constants in unlocked phase. The integration of the last 
term in a layer can be dropped to a line integral along the boundary of the layer that vanishes 
in the absence of dislocations. 

On the other hand, in locked phase the phason elastic free energy is given by 

F = \J d 3 x Ki(\V ± w 1 \ + \V ± w 2 \)+K 2 (\d z w 1 \ + \d z w 2 \) (9) 

where K\ and K 2 are elastic constants in locked phase and Vj_ = e x d x + e y d y is two 
dimensional derivative in the plane perpendicular to stacking direction (z-axis). 

To see whether a system is in a locked phase or in an unlocked phase, we could, in 
principle, compare the free energies of the systems with different uniform phason strain E 
. To do this we would need to know the energy and the entropy of the system with low 
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uniform background phason strain E, as functions of temperature. We attempted to estimate 
the entropy, using the " energy method ", "the variance method " and "the histogram 
method" but we could not estimate entropy accurately enough to distinguish the phason 
elasticity. Hence we resorted several other measures including phason fluctuations, lane 
widths and trail magnetization(defmed below). 



A. Phason Fluctuations 

The phason elastic free energy in the unlocked phase (8) can be diagonalized by Fourier 
transforming of phason variable w(x), 

w(p) = V- 1/2 J d 3 x e- ip x w(x), (10) 

where V is the volume of the system. Then free energy (8) becomes 

f = \ E km + pI)(Mp) 2 + Mp) 2 ) + K 2 pl( Wl (p) 2 + w 2 (p) 2 ). (11) 



2 

i) 

a 



IpK- 



(We ignore the third term in Eq (8) which vanishes in any configuration in our simulation. ) 
Since the free energy (11) is harmonic in Wi(p) for i = 1,2, it is straightforward to calculate 
< |w(p)| 2 >, 

<|w(P)i'>=E<Wp)i'>= JflM 4 ) + jfari . <*) 

where the angular brackets denote ensemble average. Elastic constants Ki and K 2 are 
obtained by measuring < |w(p)| 2 > for p = p x e x + p y e y and for p = p z e z . Let us define 

2 



K »<P) = TTZT-n n .M 2sn2 (13) 



< \vr(Px,P v ,0)\ 2 > (pI+pD 
2 

< |w( 0, 0,p z )\ 2 > p'j 



If the system is in the unlocked phase, lfi(p) and ^(p) should be constants, independent 
of |p| for ^ < |p| < A , where L is the system size and A is short wavelength cutoff which 
depends on the coarse graining scheme (A is order of — for a smooth function w(x), where 
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a is size of a unit cell) . A proper coarse graining scheme is necessary to get elastic constants 
by measuring < |w(p)| 2 > in our simulation. This is because the perp-space position 
(from which a smooth function w(x) is constructed) fluctuates strongly for near-neighbor 
vertices ( Ax i /Ax ~ 1 for Ax ~ 1). 

We construct the phason field w(x) from the perp-space position x- 1 by 



w x 



) = J d 3 x K(i 



x') S ± W 



(14) 



where K(x) is a smearing kernel satisfying / d 3 x K{x) = 1 and 5 ,_L (x) is the representative 
surface — piecewise linear interpolating function of xi. Precisely, we use 



^(x) = 



x i( x ) 



for x at a vertex of the upper face of a cell, 



i(^( Xl ) + ^(x 2 )) 



Xp h (x a )c4 + Xp h (xb)<i a for x in an edge (d a from one end x a , 

db from the other end x&) of the upper face, 
for x inside the upper face, where xi 
and x 2 are the projections of x to the 
edges (meet at a vertex) of the upper face, 
for x inside the unit cell, where x up is 
the projection of x to the upper face 



up J 



and 



1 1 



47r r 



e(2-r xy )5(z), 



(15) 



where r 



[x 



2 + y 2 ) 1 ^ 2 an d Q(x) is the step function. When this coarse graining scheme 



is applied to the 2D Penrose model, the known elastic constants^ are recovered for |p| < 
A ~ 1.5. Note that a simple-minded coarse graining scheme to replace Eq. (15), w(x) = 
J2iZi x p?t( x i) <H X — x «) ( w (p) = Si x p/i( x e~* p x ) would not have worked as well. In 
the 2D Penrose model, the known elastic constants^ are recovered for only |p| < A ~ 0.2 
which is order of long wavelength cutoff (^) of the biggest system in our simulation. 

To determine if the system is in the unlocked phase, we measure whether i^i(p) and 
i^2(p) a s defined in Eq. (13) are |p| - independent. In Fig. 3, we show a plot of Ki(p) 
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and K 2 (p) vs |p| at T = 1.5 and 2.5 for a sequence of lattices of increasing size. These 
plots show that ifi(p) and ^(p) are constants and independent of system size, indicating 
that the system is unlocked at T > 1.5 . We estimate the elastic constants to be K\ = 
1.72 ±0.02, K 2 = 0.56 ±0.02 at T = 1.5 and K x = 1.22 ±0.02, K 2 = 0.14 ±0.01 at T = 2.5. 

In contrast, Fig. 4 illustrates the same calculation for T = 1.0 . Here, Ki(p) and K 2 (p) 
fluctuate wildly with |p| and the mean value appears to diverge with increasing system size. 
Hence T = 1.0 is clearly below a phason transition (out of the unlocked phase). The wild 
behavior of Ki(p) and K 2 (p) is consistent with the notion that the phason elastic energy 
is locked. Hence, phase fluctuation measurements can be used to establish a transition out 
of the unlocked phase described by an elastic energy of the form in Eq. (8). However we 
cannot prove from phason fluctuation measurements that the low temperature phase has 
F ~ |Vw| (as expected for a Penrose tiling phase). This is because the phase fluctuations 
appear to become pinned (as indicated by the divergence of K\ and K 2 ). Since w is not 
thermodynamically excited, the dependence of F on |Vw| cannot be measured. 

Further evidences of the transition from the unlocked phase is provided by the measure- 
ment of the average phason field within a layer: 

w(z) = -g J s w ( x , y, z ) dxd v, ( 16 ) 

(where S is the area of a layer) and then calculating how this average fluctuates from layer 
to layer. The average phason field within a layer w(z) is related to the Fourier components 
of phason field at p x = p y = by, 

W(z) = V-^J2™(P* =Py = ^V z Y Vz \ (17) 

Pz 

since 

™(Px =P y = 0,Pz) = J w(x,y,z)e~ iPzZ dv 

= (j-f 2 I ™(z)e~ iPzZ dz, (18) 
L z J 

where L z is the number of layers. In an unlocked phase \w(p x = p y = 0,p z )\ 2 = 2/{K 2 p 2 z ) 
(Eq. (12)). Hence, the mean square fluctuation of w(z) in the unlocked phase: 
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< (aw) 2 > = -!- <EN^)-f E^)l 2 > 

" V 27T i c /L z K 2 ^ 

where c is order of 2n and c/a is upper wave number cutoff. Fig. 5 shows < (Aw) 2 > 
vs L z at T — 1, T = 1.5 and T = 2.5. The initial configuration within each layer is a 
6th approximant. < (Aw) 2 >does not show a linear dependence on L z at T = 1.0 while 
it is linear in L z at T = 1.5 and T = 2.5 (consistent with the earlier conclusion of an 
unlocked phase at T = 1.5 and T = 2.5). From the slope in Fig. 5.b and Eq. (19), we find 
cKi = 1-10 ± 0.06 at T = 1.5 and CK2 = 0.27 ± .03 at T = 2.5. Comparing these values to 
the elastic constants from the measurements of < |w(p)| 2 >, implies cutoff constant c ~ 2 
in Eq. (19) while we have a ~ 1 from Fig. 5. 



B. Lane widths 

As a mean of analyzing the low-T phase, we have measured the spacing between two 
adjacent trails in a layer. 

A "trail" in our model is a contiguous strip of tiles which share a common side face 
direction (side face direction q Q of a side face parallel to the plane spanned by and e| 
is q a = x e|, a = 0, ... ,4). We shall call a trail which has a side face direction q a , an 
a-direction trail. An a-direction trail runs along the q Q -direction on average. As shown in 
Fig. 6, each layer consists of sets of parallel trails (Fig. 6 shows a layer in a 5th approximant 
system viewed from the ten-fold axis). The regions between the trails will be referred to as 
"lanes". In a perfect Penrose tiling, two different widths of lanes exist: thick lane and thin 
lane. These lanes repeat quasiperiodically (Fibonacci sequence) in the direction normal to 
the trail direction. Hence, in a locked phase, where this quasiperiodicity is believed not to 
be destroyed, the distribution of the lane widths is bimodal. The distribution under the 
one mode (corresponding thick lanes) is r times bigger than that under the other mode 
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where r is Golden mean. In the unlocked phase, the distribution of the lane widths may 
have merged into one peak (as the distribution of the distances between nearest balls in a 
unlocked phase of Frenkel-Kontorova modeli). With this in mind, we have measured the 
number of lanes N lane P(W)dW whose average width is in between W and W + dW where 
jyiane j g total number of lanes in the system (average width of a lane is defined as the 
average of the spacing between two trails contiguous with the lane (see Fig. 6)). At low-T 
(at T = 1.0), as shown in Fig. 7. a, P(W) has two peaks near the values corresponding the 
lane widths (W[ ock and Wlf ck ) of a Penrose tiling. As the system size diverges, lane widths 
converge the values of Penrose tiling lane widths. In contrast, P(W) at high-T (at T—1.5, 
Fig. 7.c) shows a distribution with one peak at W unlock = L/N l ^ ne as L — > oo, where L is 
the system size and N l £ ne is the number of a-direction lanes in a plane (A^ ane = f k for a 
kth approximant). P(W) near the transition temperature (at T = 1.3) is shown in Fig. 7.b. 
From this figure, it is hard to tell whether the curve P(W) at T = 1.3 will be bimodal (as 
in a locked phase) or be monomodal (as in a unlocked phase) as L —>■ oo. We have checked 
reversality by doing some cooling runs from above T c . The distribution of the lane widths 
merged into one peak at high-T, come back to original bimodal mode with peaks at Penrose 
tiling lane widths. This measurement of the distributions of the lane widths, suggests that 
the system is in locked phase at low temperature. 

To estimate the transition temperature T c we have measured the number of lanes whose 
widths are around the lane widths of a Penrose tiling and lanes with widths near w unlock . 
Let us define P lock and p unlock as the relative number of these lanes: 

p iock= f P (w) dW+ f p(W)dW, 

JRi JR2 
punlock = f pf W )dW, 

p _ plock punlock (20) 



where R x = (W[ ock - 8 X , W\ ock + <$i), R 2 = (W l 2 ock - <Ji, W\ ock + 81), and R 3 = (W unlock - 
82, W unlock + 82)- Here, 81 and 82 are small numbers (<C 1) with 82 = 28\. Then P should be 
1 at T = and should be negative at T = 00 (the precise value at T = 00 depends on the 
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choice of 8\ and 82 but greater than —1 always). 

Fig. 8. shows P vs temperature for a sequence of lattices of increasing size. We have 
chosen 8\ = 0.05 and 82 = 0.1 . P converges to 1 at low T and has negative values at high T. 
We roughly estimate the transition temperature T c ~ 1.3 from Fig. 8 as the value around 
which the graphs of different sizes cross each other. 

C. Trail Magnetization 

As a new order parameter to analyze the low-T phase and the transition from the low-T 
phase, we have devised a novel measure that we have termed " trail magnetization". Trail 
magnetization traces the ordering of "hexagonal prisms" along worms and trails that run 
through each layer. 

A hexagonal prism, as shown in Fig. 2, consists of three tiles. The type (D or Q) of 
a hexagonal prism is identified with the type of the center vertex on the upper hexagon of 
the hexagonal prism0. By convention, the orientation of a hexagonal prism is "+"["—"] if 
the center vertex has an edge leading away from it along direction e| [— ejjj. A worm in a 
layer is an unbroken sequence of hexagonal prisms and the length of a worm is defined to be 
the number of consecutive, connected hexagonal prisms. Flipping one hexagonal prism in 
a perfect worm creates a mismatch along either side face adjoining the worm ( flipping the 
hexagonal prism again annihilates the mismatches). If we restrict flipping to only one worm, 
the 2D Penrose model is analogous to the one dimensional Ising model, assigning 'spin up' 
for one orientation hexagon and 'spin down' for the flipped hexagon in a worm. If the 2D 
Penrose tiling were an aggregation of uncoupled worms ( i. e uncoupled one dimensional Ising 
models), we would expect an order- disorder transition at T = which coincides with the 
result that the transition to the unlocked phase in 2D Penrose tiling is at T = 0. 

Along a trail, many worms( chains of hexagonal prisms) may be found. Along any 
trail in a perfect Penrose tiling, all worms longer than one hexagonal prism have the same 
orientation. Different trails may have opposite hexagonal prism orientation. Also, some 
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worms of length one point in the opposite direction to the other hexagons in a trail. All 
together, about 97% of hexagonal prisms in a trail have same orientation in a perfect Penrose 
tiling while the hexagonal prism in a maximally random tiling are oriented randomly. Hence, 
we define the trail magnetization in a orientation 

< = nL E I E (2i) 

a ie{tr a } jetri 

where N^ exa is number of hexagons in a-direction trails, and {tr a } means all trails (f k x L z 
trails for a kill approximant) in a-direction and tr l a means the ith trail in a-direction (i = 
1, . . . , fk x L z ). Here, spin variable Sj is assigned to the jth hexagonal prism and takes ±1 
depend on the orientation of the hexagonal prism. Without loss of generality, we will discuss 
results for a = (ejj-direction trails). 

Figure 9 shows the trail magnetization and "trail susceptibility" 

Xm = N hexa i(< (m tr ) 2 > - < m tr > 2 ), (22) 

and illustrates how m tr can serve as a useful diagnostic measure. At T = 0, the trail 
magnetization converges to a fixed value as the system diverges, a value consistent with 
the expectation value for a locked, Penrose tiling phase. At T > 1.5, the trail magneti- 
zation approaches zero as L — > oo, consistent with an unlocked phase. The magnitude of 
susceptibility maximum seems to diverge near the transition temperature. To obtain the 
transition temperature and critical exponents, we attempt to fit the trail magnetization and 
susceptibility according to 

m tr {T,L) = L-V" f[(T -T C )I}% (23) 
Xm (T,L) = L^g[(T-T c )L 1 /% (24) 

where L is the system size, T c is transition temperature of the infinite system and critical 
exponents (3, 7, v are defined from the temperature dependence of order parameter m tr , 
susceptibility Xm, and the correlation length £ near the transition temperature: 
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m* ~(T-T c )f>, for (T < T c ) 
Y ~ IT - T I" 7 

Am \ A ± c\ i 

£~\T-T C \- V . 



In Fig. lO.a, we plot m tr (T,L) L^ u vs (T - T C )L X I V and in Fig. lO.b, Xm(T,L) L~^ u vs 
(T - TJL 1 ^, for values 

T c = 1.24 /3 = 0.2 
7 = 1.6 z/ = 1.6. 

The m tr and Xm curves for various system size superimpose clearly. We have checked the 
possibility of a sequence of transitions rather than a single transition by considering 

' E I E (25) 



m 



/tr 



a AThexa 

a i€{tr a } j t etr 



where tr l a traces over the ith trails in all layers along direction a. (To see the ordering in 
stacking direction also, we took a sum of spin valuables in ith trails of all L z layers before 
taking the absolute values instead of summing the spin values in each individual trail as in 
m tr '.) We have been able to make m' tr and Xm' curves for various system size superimpose. 
The transition temperature from the measurement of m' tr is consistent with the result from 



m tr . 



We also measured energy per tile < e >=< E/N > and specific heat C v = (|;) 2 (< e 2 > 
— < e > 2 )N, where E is the system energy given by Eq. (1) and iV is number of tiles. 
Specific heat has its maximum near the transition temperature T c and the magnitude of its 
maximum seems to be independent of the system size (Figure 11) implying that the specific 
heat exponent a = where a is defined by, 

C v ~ \T-T c \- a . 

The results of trail magnetization and specific heat measurements seem very consistent 
with a single continuous transition from a locked phase {m tr 1) to an unlocked phase. 
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To check the robustness of the result, we have repeated the analysis for other related 
models. First we changed the interaction strength ratio ei/e 2 . Recall that the net energy- 
cost of a hexagonal prism flip is 2ei + 6e 2 . We have tested ei/e 2 = 1 and ei/e 2 = 9 and 
found that the systems show a locked phason to unlocked phason phase transition at finite 
temperature. 

We have also considered models where we have relaxed the interlayer "stacking con- 
straint" (that a hexagonal prism can be flipped only if a hexagonal prism lies just above 
or below with sharing a same hexagon boundary in between; see sec. II). The low temper- 
ature phase of this model appears to be the same as that of the model with the stacking 
constraint. At high temperature (T > T c ), free energy shows a quadratic dependence on 
spatial variation in phason variable in a layer (|<9 x w| 2 + l^wj 2 ). For the stacking direction, 
free energy shows a quadratic behavior in (<9 2 w) up to some temperature T' c {T' c > T c ) which 
depends on the system size. Above the T' c) layers become decoupled so that the averaged 
phason field w(z) of zth layer (Eq. (16)) is not correlated with w(z± 1) (correlation length 
of w(z) is zero for T > T' c ). However T' c increases as the system size is getting bigger and 
we speculate that T' c may diverge in the thermodynamic limits. 

V. CONCLUSIONS 

Our paper presents systems that show a finite temperature locked phason to unlocked 
phason phase transition in quasicrystals. All models for 3D decagonal quasicrystals we have 
studied show a single continuous phase transition at finite temperature from low-T locked 
phase to high-T unlocked phase. 

Phason fluctuations of the system show strong evidences that the system is in unlocked 
phase at T > T c in which free energy is described as square gradient of phason variables. 
At T < T c , from the measurements of "lane width" and "trail magnetization" we conclude 
that the system is in locked phase. The finite specific heat peak and the critical exponents 
we obtained from the scaling behaviors of the m tr show that the transition is continuous. 
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Recently Hiraga et. al. have noted that the periodically spaced layers in AlPbMn 
decagonal quasicrystals strongly correlated atomic order. Our simulations show that there 
is a single transition in which both intralayer and interlayer phason fluctuations transform 
from locked to unlocked. Hence, the observation of correlated order (locking) between layers 
would imply that that AlPbMn quasicrystals are in the locked phase and not in the random 
tiling in spite of some apparent disorder within the layer. (The disorder is likely to be due 
to decapod defects or unusual local isomorphism class.) 

To predict the transition temperature T c in real quasicrystals, we need to know how big 
the mismatch energies (ei and e 2 ) are. Note that T c is of order the geometric mean mismatch 
energy e = (eie 2 )^ for our model. T c will approach zero if either e\ or e 2 approaches zero. 

The unlocking transition could be preempted by the melting transition. If the transition 
temperature T c is higher than the melting temperature T m , quasicrystals remain in the locked 
phase for temperatures ranging all the way up to melting point. For these systems, one would 
expect more quenched phasons than for a system which has unlocked phase between the melt 
phase and the locked phase since unlocking implies rapid relaxation of phason fluctuations. 
Hence, the relation between T c and T m may partially account for the reason why some 
systems form near-perfect quasicrystals and others do not. For the system which has an 
unlocked phase (T c < T m ), the transition from the unlocked phase to the locked phase could 
be observed in experiments. An observational effect could be Debye- Waller suppression of 
the diffraction peak intensities. In the unlocked phase, phasons can be thermodynamically 
excited. As we can see in Fig. 3, the phason elastic constants increase with decreasing 
temperature in the unlocked phase. Consequently, the Debye- Waller suppression decreases 
as T decreases toward T c . Then, after the transition to the locked phase, thermal phason 
fluctuations are frozen and the the Debye- Waller suppression disappears. 
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Figure Caption 

Fig. 1 (a) Skinny Cell: Upper and lower faces are in the shape of a Penrose skinny rhom- 
bus. Side faces have arrow marks according to the arrow patterns of the Penrose skinny 
rhombus. 

(b) Fat Cell: Upper and lower faces are in the shape of a Penrose fat rhombus. Side faces 
have arrow marks according to the arrow patterns of the Penrose fat rhombus. 
Fig. 2 Hexagonal prisms: 

(a) D-type hexagonal prism before (+D) and after (-D) flip, (b) Q-type hexagonal prism 
before (+Q) and after (-Q) flip. Solid circles indicate the center vertices of the upper 
hexagons. The orientation (± sign) is given by the convention explained in the text (section 
IV. c). 

Fig. 3 The elastic constants ifi(p) (a) and ^(p) (b) defined by Eq. (13) vs the magnitude 

of the wave vector |p| with unlocked phases (at T = 1.5 and T = 2.5). ifi(p) and K 2 (p) 

are constant over |p| and independent of the system size. In this and the following figures, 

the numbers in the legend represent the numbers of tiles N in the systems. 

Fig. 4 The elastic constants ifi(p) (a) and ^(p) (b) vs the magnitude of the wave vector 

|p| with locked phase (at T = 1.0). ifi(p) and ^(p) increase in magnitude with increasing 

system size. Points of data for each system size are connected by lines shown in the legend. 

The average of K; t over |p| for each system size is indicated by an arrow. 

Fig. 5 Mean square fluctuation, (Aw) 2 , of w(z) is plotted vs stacking direction size L z 

with locked phase (a) and unlocked phases (b). 

Fig. 6 Trails and lanes in a layer (above) and average widths of lanes (below). The sequences 

of the shaded tiles are trails. The region between neighboring trails is a lane. Average lane 

width is defined as the average of the spacing between two trails contiguous with that lane. 

Fig. 7 Lane width distributions with locked phase (a) and unlocked phase (c). 

Lane width distributions near the transition temperature are shown in (b). 

Fig. 8 System size dependence of P defined in Eq. (20) vs temperature. 

Fig. 9 System size dependence of the trail magnetization (a) and susceptibility (b) plotted 
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against temperature. 

Fig. 10 Finite-size scaling plots of the data for trail magnetization, shown in Fig. 9. a, and 
for susceptibility, shown in Fig. 9.b. Here, mL 13 ^ (a) and XmL' 1 ^ (b) are plotted versus 
(T — T^L l l v with the following choice of exponents: (3 = 0.2, 7 = 1.6, and v = 1.6 and the 
transition temperature T c = 1.24. 

Fig. 11 System size dependence of specific heat plotted against temperature. 
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